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We study the Mott-insulator transition of bosonic atoms in optical lattices. Using perturbation 
theory, we analyze the deviations from the mean-field Gutzwiller ansatz, which become appreciable 
for intermediate values of the ratio between hopping amplitude and interaction energy. We discuss 
corrections to number fluctuations, order parameter, and compressibility. In particular, we improve 
the description of the short-range correlations in the one-particle density matrix. These correc- 
tions are important for experimentally observed expansion patterns, both for bulk lattices and in a 
confining trap potential. 
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I. INTRODUCTION 

The Mott-Hubbard model of interacting bosons on 
a lattice has been used to describe superfluid Mott- 
insulator transitions in a variety of systems, e.g., Joseph- 
son arrays and granular superconductors yj. The recent 
suggestion Q to experimentally observe this transition in 
a system of cold bosonic atoms in an optical lattice and 
its successful experimental demonstration Q has rekin- 
dled the interest in the Mott-insulator transition and trig- 
gered agreat deal of theoretical 0J1JE 0J1 H EH 
111 [H [H Q E3 and experimental [H 0E||i activity. 
The possibility to directly manipulate and test the many- 
body behavior of a system of trapped bosonic atoms in an 
optical lattice jlall7| is very attractive. Possible applica- 
tions include the use of a Mott state of bosonic atoms in 
an optical lattice as a starting point to create controlled 
multiparticie entanglement as an essential ingredient for 
quantum computation 0, 0, 

The Mott-insulator quantum phase transition is driven 
by the interplay of the repulsive interaction of bosons on 
the same lattice site and the kinetic energy. Hence the 
ratio of the onsite energy and the bandwidth forms the 
key parameter in the system. In optical lattices, this 
parameter can be easily controlled and varied by sev- 
eral orders of magnitude, enabling detailed studies of the 
quantum phase transition. Probing the system by taking 
absorption pictures to image the expansion patterns af- 
ter a reasonable expansion time yields information about 
the momentum distribution of the state. This procedure 
was used to experimentally confirm the Mott transition 
in an optical lattice t 3j . 

The essential physics of cold bosonic atoms in an 
optical lattice is captured by a bosonic Mott-Hubbard 



model describing the competition between hopping and 
on-site interaction. A number of approximation schemes 
have been used to study this model analytically Q, |2(], 
l2ll |2^ | as well as numerically, using approaches like 
the Gutzwiller mean-field ansatz |2j. |23| . density-matrix 
renormalization group (DMRG) |24l l25l |26| . exact di- 
agonalization ( EDl|2Ti 1281 and Quantum Monte Carlo 

(qmc) nHfllHilllila. 



In this article, we study the short-range correlations, 
not included by the Gutzwiller ansatz, by using pertur- 
bation theory. The main purpose is to find corrections 
to the short-range behavior of the one-particle density 
matrix, which is directly relevant to experimentally ob- 
served expansion patterns. These patterns are important 
for determining the location of the insulator-superfluid 
transition. We note that in the insulating state our per- 
turbative approach is identical to the one used in |2(j 
(see also |34|), although there the goal was different, viz., 
studying corrections to the phase diagram. 



The remainder of the article is organized as follows: In 
Section[Hj we will introduce the model and its mean-field 
solution. The general perturbative approach is briefly 
outlined in Section [illl while details may be found in the 
Appendix. Numerical results are presented and discussed 
in Section llVl first for local observables l|IV All and then 
for the density matrix (|IV Bfl 



Implications for expansion 
patterns both for bulk systems and a harmonic confining 
potential are discussed in Section HV CI 
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II. THE MODEL 

The cold bosonic gas in the optical lattice can be de- 
scribed by a Mott-Hubbard model Q 

M M 

h = -wni(m - 1) + Y ( e * -rfm-J^ a \ a 3 • (!) 
i=l »=1 

Here, M is the total number of lattice sites, a\ (<Zj) creates 
(annihilates) a boson on site i, = ajaj, ?7 is the on- 
site repulsion describing the interaction between bosons 
on the same lattice site, and /z denotes the chemical po- 
tential. The kinetic term includes only hopping between 
nearest-neighbor sites, this is denoted by the summation 
index J is the hopping matrix element that we 

will assume to be lattice-site independent. Finally, e$ 
describes an external on-site potential that is commonly 
present in experiments. 

The Gutzwiller (GW) approach is based on an ansatz 
for the many-body ground state that factorizes into single 
lattice-site wavefunctions 

M / oo \ 

ico>=n ■ ( 2 ) 

i=l \n=0 ) 

The Gutzwiller wavefunction represents the ground 
state of the following mean-field version of the Mott- 
Hubbard Hamiltonian, Eq. (JTJ: 

M M 

h M f = ^ 7r n *(™» - 1 ) + X! ( £i _ n% 

(3) 

Here ^ is the mean-field potential on the z-th lattice 
site, which is self-consistently defined as the expectation 
value of di in terms of the Gutzwiller wavefunction, = 

(GoHGo) m. 

Using the Gutzwiller ansatz to obtain an approximate 
variational solution for the Mott-Hubbard Hamiltonian 
||TJ corresponds, however, to restricting the Hubert space 
to the subset of product states. Consequently, even in 
higher dimensions, this ansatz fails to describe the cor- 
rect behavior of short-range correlations between differ- 
ent lattice sites, which are important for experimentally 
measurable observables, such as expansion patterns (mo- 
mentum distributions). Nevertheless, in the thermody- 
namic limit and higher dimensions, the Gutzwiller wave- 
function provides a good approximation in the limits of 
J — » and U — » (i.e., deep in the Mott insulator (MI) 
and superfluid (SF) phases). To get a satisfactory de- 
scription of the short-range correlations we will now de- 
rive perturbative corrections to the Gutzwiller mean-field 
result. 



III. PERTURBATIVE APPROACH 

Our aim is to start from the Gutzwiller approximation 
and improve it by perturbatively including the short- 
range correlations between lattice sites. We re-express 
the Mott-Hubbard Hamiltonian by adding the ap- 
propriate perturbation to the mean-field Hamiltonian, 
Eq. ©: 

H = H MF + V, (4) 

with 

V = -J^(al-^*)(a J -^ J ). (5) 

As the mean-field Hamiltonian represents a sum of single 
lattice-site Hamiltonians, the excited states \i a ) and the 
excitation spectrum e(*' Q ) can be obtained numerically 
for each lattice site i separately. Hence we can write the 
excitations of Hmf as product states of single lattice-site 
excitations, 

M 

\c a ) =niv) . (6) 

and the excitation spectrum as a sum over the single 
lattice-site excitation energies, 

M 

e a = (7) 

i=i 

where a = (cci, oli, . . .) describes the set of quantum num- 
bers characterizing the given many-body energy eigen- 
state. 

Having obtained the mean-field solution from the 
Gutzwiller ansatz, we can now proceed to improve our 
wavefunction performing Rayleigh-Schrodinger perturba- 
tion theory |37j in V: 

(G a \Gi) = 

(G"|U|G ; _ 1 )-E'rii£n(G"l^-») ( 8 ) 
eo~e a 

where |G;) and e; are the l-th order corrections to the 
wavefunction and grand-canonical energy E — fiN, re- 
spectively. 

Knowing the excited wavefunctions, Eq. JBJ, and the 
excitation spectrum, Eq. J7J), perturbative corrections to 
observables can be calculated explicitly. The resulting 
expressions up to second order in V , which we use in the 
following, are derived and illustrated diagrammatically 
in Appendix IaI 

IV. NUMERICAL RESULTS 

We start by computing the Gutzwiller wavefunction 
numerically using a conjugate-gradient descent method. 
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FIG. 1: Results from perturbation theory for homogeneous 
lattices in d = 1 (diamonds), d — 2 (triangles), and d = 3 
(squares) dimensions. The dashed-dotted lines in Fig. 
c are the Gutzwiller results, (a) Number fluctuations <7j = 
V ( n i ) ( 71-i ) 2 calculated for a commensurate filling of one 
boson per lattice site. The dashed line shows the result from 
the exact diagonalization for 7 lattice sites and N — 7 bosons, 
(b) Order parameter $i = (en) and (c) compressibility k both 
computed at a fixed chemical potential fi/U = 0.5. 



Propagation steps of the time-dependent Gutzwiller 
equations 2] in imaginary time are performed to test 
whether the minimum found by the conjugate-gradient 
descent method is indeed the ground state. 

Afterwards, we calculate eigenenergies and -vectors of 
each site, with respect to the mean-field Hamiltonian. 
As explained above, this forms the basis for perturba- 
tion theory, which we use to calculate the corrections to 
the observables up to second order. The results obtained 
from perturbation theory show important modifications 
to single-site observables as well as to the correlation 
function. 



A. Single Lattice-Site Observables 

These observables are composed of operators acting 
only on single lattice sites. Hence, they describe local 
properties and will be less sensitive to the correlations 
between lattice sites. Thus, values for single-site observ- 
ables obtained from the Gutzwiller wavefunction already 
provide a good approximation in most cases (unless fluc- 
tuations are concerned). 

The leading corrections to the mean values of single 
lattice-site observables (SLSO) are of second order in V. 
The results for the order parameter $i , the compressibil- 
ity k and the number fluctuations are shown in Fig. ^ 
where the ration J/U has been scaled by the dimension, 
to keep the same mean- field transition point. The solid 
lines in Fig. ^ show the results from perturbation the- 
ory, as compared to the Gutzwiller result (shown as the 
dashed - dotted line) . All three quantities show a vanish- 
ing perturbative correction both for small and large J/U, 



as the Gutzwiller wavefunction becomes a good approx- 
imation in these regimes (for lattice dimensions d > 1). 
As expected, deviations from the mean-field picture are 
strongest near the MI-SF transition, where higher-order 
corrections will become more and more important. 

The order parameter $i — (aj) shown in Fig. \T]p 
gets suppressed in the SF. Perturbative corrections to 
the Gutzwiller result are particularly large in ID (where 
the order parameter vanishes in reality), but get smaller 
with increasing dimension. This is not surprising, as 
Gutzwiller is a mean-field approach and hence a bet- 
ter approximation for higher-dimensional systems. The 
critical value (J/U) c is not modified within the present 
perturbative approach. 

Figure shows the results for the compressibility 



M ^ d{m) 

2 2^ 



N 



dfi 



(9) 



The results of perturbation theory (PT) show a decrease 
of the compressibility, pointing to an increasing stiffness 
of the SF phase induced by the short-range interaction. 

Finally, we computed the local particle number fluctua- 
tions Oj = \J (nf) — (rii) 2 . Exact diagonalization calcula- 
tions have shown that the number fluctuations are chang- 
ing smoothly at the MI-SF transition [^3, whereas the 
Gutzwiller result predicts vanishing fluctuations in the 
MI phase, cr^ = (see dashed-dotted line in Fig. HJ;) . Our 
perturbative results reproduce the non-vanishing part of 
CTj in the MI regime and agree well with our result ob- 
tained from exact diagonalization of a one-dimensional 
system of 7 lattice sites. Significant deviations from the 
exact ID result are seen in the MI-SF transition regions, 
starting from the mean-field critical value, (J/U) c = 
0.086, up to values of the order of the critical values of 
(J/U) c w 0.2 ■•■0.3 usually obtained from DMRG HI 
and QMC [3(| calculations. 

Nevertheless, based on the good agreement of our per- 
turbative result with the exact diagonalization for a large 
region in the MI, we conclude that the number fluctua- 
tions are mainly produced by next-neighbor particle-hole 
fluctuations included in perturbation theory. 



B. The Correlation Function 

The single-particle density matrix py L = (a\a,j) is of 
particular importance as it describes the correlation be- 
tween the different lattice sites. The correlation function 
Pij shows off-diagonal long-range order in the SF state 
(in dimensions d > 1), in contrast to the MI phase where 
Pij decays exponentially. 

The experimental observation Q of the MI transition 
relies on the different behavior of the density matrix in 
the MI and SF regimes, which can be visualized by taking 
absorption pictures of the freely expanding atomic cloud. 
Assuming that the expansion time is long enough and 
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FIG. 2: Short-range behavior of the one-particle density ma- 
trix pij as a function of site distance \i — j\. Results for a 
homogeneous lattice of 7 sites with N = 7 Bosons have been 
obtained from exact diagonalization (ED), from second-order 
perturbation theory (PT) , and from the Gutzwiller mean-field 
ansatz (GW). (a) J/U = 0.05 (MI regime), (b) J/U = 0.1, 
(c) J/U = 0.2, (d) J/U = 0.5 (deep SF regime). All calcula- 
tions use periodic boundary conditions. The mean-field value 
(J/U) c for the MI transition in the commensurate case with 
one boson per lattice site is (J/U) c = 0.086 in ID (see Fig.0. 
(The MF value differs strongly from (J/U) c = 0.277 derived 
from QMC calculations or {J/U) c = 0.260 for DMRG 
calculations PH - ). 



that the gas is dilute enough (such that atom-atom inter- 
actions can be neglected during the expansion), the shape 
of the cloud reflects the initial momentum distribution 
p(k), which is directly given by the Fourier-transform of 
the density matrix p t j : 



P (k) = Hk) 



M 

X! P"' 



(10) 



Here, w(k) is the Fourier transform of the Wannier 
functions w(r — r^) describing the wavefunction of a sin- 
gle lattice site. The presence of the factor w(k) in 
Eq. (|10|) provides a cutoff at high momenta. 

The mean-field results for py only describe the differ- 
ent long-range behavior in the MI and SF. For a homoge- 
neous lattice the correlation function calculated from the 
Gutzwiller wavefunction gives pa = rii for the diagonal 
elements and then drops instantly to = |<I>| 2 for all 
off-diagonal elements i ^ j. Short-range correlations are 
not reproduced by the Gutzwiller approach. This devia- 
tion is particularly severe in the MI, where the mean-field 
result predicts a completely flat momentum distribution, 
whereas the short-range correlations (i.e. the exponen- 
tial decay of py) yield smooth bumps in the expansion 
pattern. These can be distinguished from the (S-peaks of 
the SF only after a sufficiently long expansion time. 



Applying perturbation theory to the GW wavefunction 
improves the structureless GW correlation function. In 
Fig. [3 we have compared the results of GW mean-field 
theory, of PT, and of exact diagonalization (ED). The di- 
agonalization has been carried out for a small ID lattice, 
where it is easily feasible. Although there is no long- 
range order in the SF phase for the ID case, where the 
density matrix exhibits a power-law decay towards zero, 
it is still reasonable to compare the short-range correla- 
tions. Indeed, we find a nice agreement between pertur- 
bation theory and exact results, not only for the MI (see 
Fig- [20, but also for the short-range behavior in the SF. 
This agreement is made possible by the fact that p^ de- 
cays only slowly and higher-order corrections show only 
negligible corrections for small lattices. An example for 
the SF case is shown in Fig.[2i. However, there is still a 
considerable difference for in Fig. [5Jl. This is not 

surprising as we do PT up to second order. Hence cor- 
relations over a distance of three and more lattice sites 
are only corrected by the global mean-field correction for 
the infinite lattice (see Eq. I[A11J) . Eq. (|A12j) . and Fig. ED 
in Appendix 0. We expect better agreement for larger 
lattice sizes as finite-size effects, arising in small lattices, 
are still considerable for M = 7 sites used in our ED 
calculations. 

Finally, for intermediate values of J/U, shown in 
Fig. Gh>,c, we observe a faster drop in the off-diagonal 
correlations, such that higher-order contributions in the 
PT become more important. In any dimension, the per- 
turbation V is no longer small at the tip of the MI-SF 
transition lobe (Fig. |2t), and the PT breaks down. Nev- 
ertheless, comparing with exact diagonalization results, 
Fig. 13) shows still good agreement with ED, in contrast 
to Fig.[2t, which shows clear disagreement. Even though 
the parameters J/U = 0.2 chosen for Fig. 0: are close 
to the MI-SF transition for ID-lattices (as predicted by 
DMRG and QMC calculations), PT reproduces 
the correct slope for the off-diagonal decay and lacks only 
the wrong offset from the mean field. Thus, even for 
this case, PT represents a qualitative improvement on 
the Gutzwiller result. The results for both approxima- 
tions, Gutzwiller and PT, are expected to become better 
in higher dimensions (with the perturbative corrections 
diminishing in size). 

As discussed before, the perturbative enhancement of 
the description of short-range correlations is expected to 
lead to strong consequences for the momentum distri- 
butions. As an example we discuss a set of momentum 
distributions p(k)/|w(k)| 2 for a homogeneous 2D lattice, 
Fig. int-fi and compare the Gutzwiller results to those 
improved by PT. The improved PT versions, Figs.|3Ji-f, 
show much finer structures than the mean-field results, 
Fig. EK-C- PT predicts broad peaks in the MI regions 
down to very small values of J/U, Fig. ®, whereas the 
Gutzwiller result without PT shows a structureless flat 
distribution for the whole MI region, Fig.|^i,b. Naturally, 
the modifications of p(k)/|w(k)| 2 are strongest near the 
phase transition, Figs.|3t and|3i. 
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FIG. 3: The central figure shows the correction to pi,i+i in 
second-order PT as a function of /j,/U and J/U. The or- 
der parameter (<f>) vanishes inside the Mott-insulating lobes, 
whose mean-field boundaries are given by the white dashed 
line. Plots (a) - (f) display the resulting momentum distri- 
bution without the Wannier form factor, p(k)/|iu(k)| 2 , calcu- 
lated for a 2D lattice with 25 x 25 lattice sites, (a)-(c) are the 
Gutzwiller mean-field results, and (d)-(f) are calculated using 
second-order PT. The inset in (c),(f) shows a cut of one peak 
taken along k y a = 0; dashed line for GW and solid line for 
the PT result. Arrows indicate the position of the respective 
plots in the (fi/U, J/U) phase diagram. The parameters used 
are: fi/U = 1.5, J/U = 0.0225 for (a) and (d); fi/U = 0.75, 
J/U = 0.01 for (b) and (e); and n/U = 0.5, J/U = 0.044 for 
(c) and (f). The gray-scales of plots belonging to the same pa- 
rameter set are identical. Expansion patterns (a),(b),(d),and 
(e) are normalized to the peak maximum; (c) and (f) are nor- 
malized to 1/20 of the peak maximum. 




FIG. 4: Occupation number m on a plane through the trap 
center for a 3D lattice with 15 3 lattice sites, (a) fi/U = 0.4, 
J/U = 0.005, and a/U = 0.02. (b) fx/U = 1.5, J/U = 0.0075, 
and a/U = 0.03. (c) fi/U = 0.4, J/U = 0.0075, and a/U = 
0.02. 



Going towards larger values J/U into the SF phase, PT 
gives rise to a suppression of the peaks (inset in Fig.[3Ji,f). 
This suppression can be larger than 20% of the original 
peak height and stems from the corrections to the mean 
field. Additionally, Fig. |2F shows broad peaks induced by 
the inclusion of short-range correlations. However, for 
large lattices, these broad peaks are small compared to 
the (finite-size broadened) SF <5-peaks. 



C. Harmonic Trap Potential 

In contrast to what was assumed in the last section, 
optical lattices used in experiments are not homogeneous. 
Magnetic or optical trapping potentials are used 0, |3(J 
to confine the atomic gas to a finite volume. 

The inhomogeneity caused by the trapping potential 
leads to slowly varying on-site energies, e*, in the Mott- 
Hubbard model Eq. (|TJ, that can be interpreted as a 
spatially varying chemical potential ^i OC ai = M — £j ■ Con- 
sequently the lattice is in general not in a pure MI or SF 
phase, but shows alternating shells of SF and MI regions. 
An example for a SF region surrounded by a MI shell is 
shown in Fig. 0p. 

Considering the slowly varying on-site energy as a spa- 
tially varying chemical potential gives a qualitative un- 
derstanding of the shell structures. The spatial variation 
of the chemical potential corresponds to a path parallel 
to the n/U axis in the /u/ZJ-J/JJ-diagram. Starting with 
the potential minimum, in the trap center, and then mov- 
ing off the center, decreases the effective local chemical 
potential. Whenever a MI-SF (SF-MI) phase boundary 
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ou/U=0.4, J/U=0.005, a/U = 0.02 
□ |i/U=1.5, J/U=0.0075, a/U=0.03 
HAJ=0.4, J/U=0.0075, a/U=0.02 




FIG. 5: Expansion pattern for a 3D-lattice with 15 -sites in 
the presence of a harmonic potential, (a) Momentum distri- 
bution without the Wannier form factor, p(k)/jui(k) | 2 , cal- 
culated along the k x direction. Filled symbols are the PT 
results and open symbols are the GW results. The graphs for 
the MI surrounded by a SF shell (diamonds) are rescaled by 
a factor 1/25. The results denoted by circles correspond to 
the occupation distribution of Fig. diamonds to Fig. |1J;, 
and squares to Fig. IUd. (b) Difference between PT and GW 
results. 



is hit along the path in the fi/U- J/[/-diagram, a change 
from a MI to a SF (SF to MI) shell appears. 

The inclusion of short-range correlations gives rise to 
considerable modifications also in the presence of an in- 
homogeneous trapping potential. Examples, calculated 
for a 3D-lattice, are shown in Fig. EJ where an underly- 
ing harmonic potential 



0=1 



(ii) 



was chosen. 

The three situations considered here correspond to a 
case with a large MI fraction (Fig.^Ji), a SF island sur- 
rounded by a MI shell (Fig. 0Jd), and a MI island sur- 
rounded by a SF phase (Fig. 0J:). Calculating the ex- 
pansion patterns for these situations, Fig. shows that 
the perturbative corrections arising from the short-range 
correlations lead to substantially different behavior in the 



different cases. 

For the almost complete MI state we get a correction 
to all wavevectors fc, with a fast drop at values close to 
the peak center k = 0, which leads to a peak broadening 
in the expansion picture (see circles in Fig. . 

Particularly large changes were found for the case of 
a SF island surrounded by a MI phase, Fig. (squares) . 
Again, corrections arise for all wavevectors, but, in con- 
trast to the almost homogeneous case, the largest in- 
crease is now found for k = 0, with changes of about 
20% of the peak-maximum. 

Finally, the reversed situation, a MI island surrounded 
by a SF phase (diamonds in Fig. does not show an 
increase of its maximum peak height but a considerable 
reduction (over 5% of the peak maximum). This is not 
surprising, as the majority of the lattice sites are now 
contributing to the SF phase, and a peak reduction was 
also observed for the bulk SF phase. 

We would now like to compare our approach with the 
results obtained by Kashurnikov et al. 33], who used 
QMC calculations to calculate expansion patterns for a 
small 3D lattice with harmonic confinement. There is 
good qualitative agreement in all cases with high su- 
perfluid fraction (compare for example Fig. |B^(b) and 
Fig. 2b(c) in Ref. 33]). However, the features of these 
expansion patterns are already well reproduced using 
the Gutzwillcr mean-field ansatz alone, in particular, the 
satellite peak, which was discussed as a signature of the 
MI-SF shell structure in [2*3 . Corrections arising from 
PT show a suppression of the SF peak, as was discussed 
above. 

However, there are also considerable discrepancies to 
the QMC results for situations with a large MI fraction, 
even after implementing second-order PT. In these cases 
(Fig.Et(d) and Fig. 2d(e) in Ref. 0), the influence of 
the MI phase on the expansion picture broadens the peak 
and leads to a homogeneous background. The discrep- 
ancies to the QMC are clearly visible in Fig. [U showing 
no peak broadening and a sattelite peak in contrast to 
the QMC results (Fig. 2d in pf). Including the short 
range correlations perturbatively corrects the expansion 
pattern in the right direction, giving rise to a suppression 
of the SF-peak. Considering the expansion pattern with 
the clearest MI features (FigUl and Fig. 2e in Ref. pf). 
we obtain the correct peak broadening from GW/PT cal- 
culation, but a larger ratio of MI background to SF peak. 

Concerning these discrepancies, we note that the ex- 
pansion pattern is highly sensitive to the value of the 
mean field. Even small deviations can lead to a change 
in the SF-peak height sufficient to mask the flat distri- 
bution of the MI phase. We checked, however, that the 
observed discrepancy is not due to a lack in accuracy 
of our numerical calculations. We therefore believe that 
the discrepancies between QMC and GW/PT for situa- 
tions with a large MI fraction can be attributed to the 
insufficiency of GW and low-order PT in describing the 
long-range correlations in this inhomogeneous situation. 
In addition, we note that the lattice employed in pi ] is 
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FIG. 6: Gutzwiller (dashed line) and perturbation theory 
(solid line) expansion patterns (cuts along k x ) for a system 
of 16 3 lattice sites considered in 0. All plots are with- 
out the Wannier form factor, p(k)/|u>(k)| 2 . The insets show 
the occupation number m for a cut along the indirection, 
(a) J/U = 1/32, (j, = 0.3775, a = 0.00610. (b) J/U = 1/80, 
H = 0.3125, a = 0.01221. (c) J/U = 1/80, fj, = 0.625, 
a = 0.01288. (d) J/U = 1/80, n = La = 0.02505. (Note 
however, that the Hamiltonian used in ; 33j differs from Eq. Q 
and hence the parameters are converted to the corresponding 
quantities used in our definitions.) 



comparatively small for the given harmonic confinement 
potential (with no complete shell of empty sites at the 
perimeter, see insets of Fig.|SJ), i.e., the choice of bound- 
ary conditions (periodic in the case of our numerical cal- 
culations) may have non-negligible effects on the outer 
lattice sites. 
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FIG. 7: Schematic diagram illustrating the term (c}cj^Vij)o 
appearing in the first-order correction (c\cj) l of the density 
matrix, where Vij = ^(cjcy + c^a) is the perturbation con- 
necting sites i and j. The diagram shows the lattice sites in 
the horizontal direction. The different steps needed to obtain 
the matrix element are shown vertically. Open circles denote 
the ground state (GS) of the given lattice site, while filled 
circles are excited states (ES) of this site, for the mean-field 
Hamiltonian Hmf- 



corrections to the correlation function, and thus to the 
expansion patterns. We compared the correlation func- 
tion obtained from PT with calculations from ED for 
small ID-lattices and found good agreement. Studying 
the expansion patterns showed that the inclusion of the 
short-range correlations to the MF-ansatz gives rise to 
distinct modifications. A broad peak can be seen in the 
PT results for the MI regime, visible even down to small 
values of J/U. Comparing PT and MF expansion pat- 
terns obtained for parameters in the SF region displayed 
a considerable suppression of the SF peak in the PT re- 
sults. Additionally, on approaching the SF-MI transi- 
tion from the SF side, broad peaks underlying the SF 
peaks were found in the PT-expansion patterns. Includ- 
ing a harmonic confinement potential leads to situations 
where SF and MI regions coexist. Hence the perturba- 
tive corrections to the expansion pattern become more 
complex. We studied lattices with different underlying 
harmonic traps giving rise to different constellations of 
SF-MI regions, finding modifications of up to 20% of the 
peak maximum in the expansion patterns. 



V. CONCLUSIONS 

We have used perturbation theory to incorporate the 
effects of short-range correlations on top of the mean 
field solution of the bosonic Mott Hubbard model. We 
derived corrections to local quantities, as well as to ob- 
servable expansion patterns. We numerically calculated 
the corrections to the MF-result, using PT up to sec- 
ond order, thus including correlations between next and 
next-nearest lattice sites. Modifications to the particle 
number fluctuations <7j, arising from the PT, gave rise to 
the expected smooth transition of Oi at the MI-SF tran- 
sition. Moreover, comparing the PT results to the ED 
results for <7j in ID lattices showed good agreement for 
small values of J/U. Of particular importance are the 
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APPENDIX A: DERIVATION OF THE MATRIX 
ELEMENTS 

We use standard stationary perturbation theory to 
calculate the corrections to the mean-field results induced 
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by the perturbation 



J 



v = -±^Vij = -jLfe + C I C ») > ( A1 ) 

where we introduce the new operators 

Ci =a,i- (ai)o = a t - ^ , (A2) 

The expectation value (•)□ is taken with respect to the 
MF wavefunction \Go). Defining the operator for the 
energy denominator as 



A eo-Hi 



MF 



(A3) 



the expectation value of an observable (A) including all 
corrections up to second order is 



(A) *(A) + (V±A)o + (A^V)o 

+ <yiv±A) + <y±A±v) 



(A4) 



1 „ 1 



1 



+ ( A & V A V) °- {V ^ V)a{A)o - 



The first line in Eq. (£Q is the MF-result (A) , followed 
by two contributions which are the first-order corrections. 
Lines two and three in Eq. (|A4(1 are the second-order 
corrections to the mean value. 



1. Density Matrix: First-Order Corrections 

As an example we will discuss the corrections to the 
density matrix = (a\a,j). The first-order correction 
to the density matrix is: 

= (cfah + + (4)1*, + ■ (as) 

For two different lattice sites i ^ j, we find 

( c l c i^ c i)o = (cj)o(4-^Ci) = , (A6) 

since the Gutzwiller ground state is a product state. As 
a consequence, the contributions (cj)i and (cj-)i vanish. 
The only remaining contributions to (ojfflj)i stem from 

(4 c i)i = ~(Vij^4 c j)o - (4 c 3-£ v v)o , (A7) 

and 

To get a better idea of the character of the terms 
arising in the perturbative expansion, we introduce, in 
Fig. \7\ a graphical representation (for the example of a 
ID-lattice). 

The graph shows a decomposition of the matrix ele- 
ment, with each row showing the wavefunction at an in- 
termediate step in the evaluation of the matrix element. 
As we deal with product states 



M 



iG«)=ni^>> 



(A8) 



we represent the wavefunction by a row of circles, where 
each circle denotes the state \i ai ) of a particular lattice 
site i. 

Open circles in Fig.[7]denote a lattice site in its ground 
state (GS), filled circles refer to an excited state (ES) 
of this particular lattice site, with respect to the local 
mean-field Hamiltonian. Note that, in general, this can 
be an arbitrarily highly excited state (although higher 
contributions are suppressed by the energy denominator, 
and a cutoff is used in practice). 

Starting with a row of open circles, denoting the GS, 
| Go), each following row corresponds to the state after 
the action of V or c\cj, as indicated on the left side of 
the graph. As all matrix elements in Eq. 1A4I) can be 
expressed in terms of a GS expectation value (-)o and a 
sequence of V and c\cj operators, the first and last row 
must always be a line of open circles. 

Let us consider for instance the second term in 
Eq. llA"7l) 



(A9) 



Reading the graph in Fig. \7\ from top to bottom corre- 
sponds to reading the matrix element from right to left. 
Starting with the GS, \Gq), the first row consists of open 
circles. The second row shows the state after the action 
of the perturbation V. Acting with Vij to the right onto 
the GS results in a state 



VijlGo) 



fa,p\ia,jp) 



(A10) 



where i and j are neighboring lattice sites and \i a ,jfj) 
denotes the state with lattice site i (j) in the exited state 
a ((3) and all other sites in their GS. Thus the second row 
shows the lattice sites i and i+ 1 in an excited state (filled 
circle), as the perturbation, Eq. IjAljl . allows only next 
neighbor interactions. Finally, the action of c\cj has to 
bring the excited states back to the GS, in order to get a 
non-vanishing contribution. Therefore, in first order PT, 
only next neighbor corrections to the correlation function 
arise, as the final row must represent the ground state 
(Go | again. The graph representing the remaining first 
term in Eq. (|A7|) is obtained by rotating the graph in 
Fig.Qby^. 



2. Density Matrix: Second-Order Corrections 

Rewriting the second-order corrections to the density 
matrix, (a j 0^)2, in terms of the operators c] and Cj gives 
Eq. i|A5() . but with (-)i replaced by (^a- In contrast to 
the first-order corrections, now the terms proportional 
to (cj)2 and (c|)2 also give non- vanishing contributions. 
Using Eq. QA4|) we obtain: 



in 



i=l 



EL(^*s^*s<V>o +T,' k {cj J kVjkiVik}o (All) 
+ E'k(V J kic J ±V jk ) . (A12) 
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FIG. 8: Diagrammatic representation of the contributions 
from \P*(cj)a and (cj)2^j to the second-order correction of 
the density matrix, (a) Terms representing contributions of 
the form Eq. 1A1H . (b) contributions of the form Eq. 1A121 . 
All notations are the same as in Fig. 
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FIG. 9: Graphs showing second-order corrections arising 
from (cjcj)2. Diagrams (a) and (b) are coming from direct- 
neighbor contributions as given by Eq. iA"l3l and Eq. llATIt 
respectively. Diagrams (c) and (d) are next-nearest-neighbor 
contributions: (c) corresponds to Eq. ljA16^ and Eq. I|A170 : 
(d) corresponds to Eq. 1A18> . 



Here, the primed sums run over all neighbors k to site 
j. The corresponding subset of graphs for ^*(cj)2 and 



{c\)2 ^j ar e given by Fig.|SJi and Fig. 03 for Eq. I|A11|I and 
Eq. (| A 12(1 respectively. Note that all terms of Eq. (|A11() 
and Eq. I|A12(I give a correction to all matrix elements of 
the density matrix independent of the distance between 
the lattice sites. We can understand these contributions 
as a modification to the MF value of the density matrix. 

For the second-order contribution, (cjc,-)2, we have to 
distinguish two cases: 

(a) Lattice site i and j being direct neighbors. In this 
case we get 

(chh = (Vid V ijhh)o + {cUd V id V ij)o (A13) 
+ <^S c ^s^)o- (A14) 



Corrections for Eq. i|A13(l and Eq. ilA14(l are shown 
in Fig. [2K an d Fig. EJd, respectively. 

(b) Configurations corresponding to two lattice sites 
i =/= j connected by two successive hopping steps 
via site k. This gives rise to six contributions: 

(chh = (A15) 

£<v; fc K Vfci i c ' Cj)o + E^4 1/jfc z c - Cj)o (A16) 

c c 

c c 
+ £<Vj*icJc, iv w )o + £ (V*i\ c h\v jk ) . (A18) 



An example for the contributions arising from 
Eq. l|AT6)l and Eq. ||aT7|) is shown in Fig.GJ. The 

last term, Eq. I|A18(1 , has the representation shown 
in Fig. [Ul. Note that for lattices with dimensions 
d > 1, the sites i,j and k need not necessarily form 
a straight line but can form a chevron. 
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